/*
  Copyright (C) 2011, Kyungjoo Kim
  
  This file is part of LINAL (LINear ALgebra)
  
  All rights reserved.

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.
  
  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.
  
  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  Also add information on how to contact you by electronic and paper mail.
  
  
  Kyungjoo Kim
  iamkyungjoo@gmail.com
*/

#include "environment.hxx"

int main(int argc, char **argv) {

  linal::initialize( argc, argv );

  bool help;
  int nthreads, datatype, m, n, b_mn;
  int
    locker_host, b_mn_host, balancer_host,
    locker_gpu, b_mn_gpu, balancer_gpu,
    n_gpu, cache, n_work, n_bin,
    com;

  double tol;

  linal::get_option_bool  ( "-help", "list options",
                            __FILE__,LINAL_OPTION_MODE_RUN, help);

  linal::get_option_int   ( "-datatype", LINAL_OPTION_TEXT_DATATYPE,
                            __FILE__,TEST_DATATYPE, datatype, help);

  linal::get_option_int   ( "-M", "size of matrix",
                            __FILE__,N, m, help);
  linal::get_option_int   ( "-N", "size of matrix",
                            __FILE__,N, n, help);
  linal::get_option_int   ( "-b_mn", "block size of matrix",
                            __FILE__,B_MN, b_mn, help);

  linal::get_option_int   ( "-nthreads", "# of threads",
                            __FILE__,1, nthreads, help);

  linal::get_option_double( "-tol", "error tolerence",
                            __FILE__,LINAL_ERROR_TOL, tol, help);


  linal::get_option_int   ( "-locker_host", "locker size",
                            __FILE__,12, locker_host, help);
  linal::get_option_int   ( "-balancer_host", "balancer ",
                            __FILE__,6, balancer_host, help);
  linal::get_option_int   ("-b_mn_host", "b_mn for host",
                           __FILE__,b_mn, b_mn_host, help);

  linal::get_option_int   ("-n_gpu", "# of gpus",
                           __FILE__,1, n_gpu, help);
  linal::get_option_int   ( "-locker_gpu", "locker size",
                            __FILE__,12, locker_gpu, help);
  linal::get_option_int   ( "-balancer_gpu", "balancer ",
                            __FILE__,6, balancer_gpu, help);
  linal::get_option_int   ("-b_mn_gpu", "b_mn for gpu",
                           __FILE__,b_mn, b_mn_gpu, help);

  linal::get_option_int   ("-cache", "cache policy for gpus: 0) no use, 1) read only \
2) full",
                           __FILE__,2, cache, help);
  linal::get_option_int   ("-n_work", "# of work cache",
                           __FILE__,3, n_work, help);
  linal::get_option_int   ("-n_bin", "# of bin",
                           __FILE__,80, n_bin, help);

  linal::get_option_int   ("-computing", "computing model",
                           __FILE__,LINAL_COM_DEVICE, com, help);

  if (help == LINAL_OPTION_MODE_DRY) {
    linal::finalize();
    return 0;
  }

  omp_set_num_threads(nthreads);

  // host device set up
  linal::host::Manager host = new linal::host::Manager_();  // group 0
  host->initialize( balancer_host );

  // gpu device set up
  std::vector<linal::gpu::Manager> gpus;                    // group 1
  for (int i=0;i<n_gpu;++i) {
    linal::gpu::Manager gpu = new linal::gpu::Manager_(i,1);
    gpu->initialize( balancer_gpu, cache, n_bin, n_work );
    gpu->create_obj( datatype, b_mn_gpu, b_mn_gpu );
    gpus.push_back( gpu );
  }

  // dispatcher set up
  linal::Dispatcher dispatcher = new linal::Dispatcher_();
  dispatcher->add_device( host );
  for (int i=0;i<n_gpu;++i)
    dispatcher->add_device( gpus[i] );


  linal::Dispatcher_::Group_ grp_host(b_mn_host, locker_host);
  dispatcher->add_group( grp_host );

  if (n_gpu) {
    linal::Dispatcher_::Group_ grp_gpu(b_mn_gpu, locker_gpu);
    dispatcher->add_group( grp_gpu );
  }

  dispatcher->lock();

  linal::set_dispatcher( dispatcher );

  std::cout << *dispatcher << std::endl;

  linal::Flat_ A[2], norm;
  linal::Hier_ HA;

  for (int i=0;i<2;++i) 
    A[i].create( datatype, m, n );

  HA.create( A[0], b_mn, b_mn );
  norm.create( linal::get_datatype_real(datatype), 1, 1 );

  int uplo = FLA_LOWER_TRIANGULAR;
  FLA_Random_spd_matrix( uplo, ~A[0] );
  FLA_Hermitianize( uplo, ~A[0] );
  FLA_Copy( ~A[0], ~A[1] );

  if (linal::get_verbose()) {
    FLA_Obj_show( "= A[0] =", ~A[0], "%4.3f", "===" );
    FLA_Obj_show( "= A[1] =", ~A[1], "%4.3f", "===" );
  }

  double t_base, t_linal, t_fla;
  double flop = linal::get_flop_lu( 0, m, n );

  t_base = FLA_Clock();
  FLA_LU_nopiv( ~A[1] );
  t_fla = FLA_Clock() - t_base;

  printf("FLA DONE\n");

  dispatcher->set_computing_model(com);
  dispatcher->stat_begin();

  t_base = FLA_Clock();
#pragma omp parallel
  {
#pragma omp single nowait
    linal::lu_nopiv( HA );
  }
  t_linal = FLA_Clock() - t_base;

  dispatcher->stat_end();
  dispatcher->set_computing_model(LINAL_COM_HOST);

  dispatcher->statistics(std::cout);

  printf("Time LINAL % E , FLA % E\n", t_linal, t_fla);
  printf("FLOP LINAL % E , FLA % E\n", flop/t_linal, flop/t_fla);


  if (linal::get_verbose()) {
    FLA_Obj_show( "= A[0] =", ~A[0], "%4.3f", "===" );
    FLA_Obj_show( "= A[1] =", ~A[1], "%4.3f", "===" );
  }

  FLA_Axpy( FLA_MINUS_ONE, ~A[0], ~A[1] );
  FLA_Norm1( ~A[1], ~norm );

  if (linal::get_verbose())
    FLA_Obj_show( "=Diff A[1]=", ~A[1], "%4.3f", "===" );

  float  norm_f = *((float *)FLA_Obj_buffer_at_view( ~norm ));
  double norm_d = *((double*)FLA_Obj_buffer_at_view( ~norm ));

  printf("- TEST::");
  for (int i=0;i<argc;++i)
    printf(" %s ", argv[i]);
  printf("\n");

  int rval;
  switch (linal::get_datatype_real(datatype)) {
  case LINAL_SINGLE_REAL:
    printf("SINGLE PRECISION\n");
    if (norm_f < LINAL_ERROR_TOL) {
      printf("PASS::Norm :: %10.9f \n", norm_f);   rval = 0;
    } else {
      printf("FAIL::Norm :: %10.9f \n", norm_f);   rval = -1;
    }
    break;
  case LINAL_DOUBLE_REAL:
    printf("DOUBLE PRECISION\n");
    if (norm_d < LINAL_ERROR_TOL) {
      printf("PASS::Norm :: %E \n", norm_d);   rval = 0;
    } else {
      printf("FAIL::Norm :: %E \n", norm_d);   rval = -1;
    }
    break;
  }

  // ---------------------------------------
  // ** Matrices
  HA.free();

  for (int i=0;i<2;++i) 
    A[i].free();
  
  norm.free();

  for (int i=0;i<n_gpu;++i) {
    linal::gpu::Manager gpu = gpus[i];
    gpu->free_obj();
    delete gpu;

  }
  delete host;

  delete dispatcher;

  linal::finalize();

  return rval;
}
